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Abstract. We report recent progress on causal viscous hydrodynamics for rel- 
ativistic heavy ion collisions. For fixed specific shear viscosity 77/s, uncertainties 
in the elliptic flow arising from initial conditions, equation of state, bulk viscosity 
and numerical viscosity, and the treatment of the highly viscous hadronic stage 
and freeze-out procedure are analysed. A comparison of current viscous hydrody- 
namic results with experimental data yields a robust upper limit 77/s < 5 X i. 

1. Introduction 

The viscosity of the quark-gluon plasma is presently a hotly debated subject. Its 
computation from first principles is difficult. It is thus desirable to try extracting it 
from experimental data. Viscous hydrodynamics provides a tool that can be used to 
attack this problem while simultaneously extending the region of applicability of the 
hydrodynamic approach beyond that of ideal fluid dynamics. 

During the last years, several groups have independently developed (2+1)- 
dimcnsional viscous hydrodynamic codes to describe the mid-rapidity spectra from 
heavy-ion collisions at Relativistic Heavy Ion Collider (RHIC) energies. First results 
were published in Refs. [1, 2, 3, 4, 5, 6, 7, 8, 9]. It was found that shear viscosity 
decelerates the longitudinal expansion, thereby initially slowing down the cooling 
process, but accelerates the transverse expansion, resulting in more radial flow and 
flatter hadron px-spectra. In non-central collisions, the elliptic flow coefficient vi was 
found to be very sensitive to shear viscosity: given the large expansion rates of heavy- 
ion collision fireballs, even minimal viscosity saturating the KSS bound 77/s > 1/Att 
[10] for the viscosity to entropy density ratio can lead to a strong (and thus easily 
measurable) suppression of vi. Assuming the availability of a well-established ideal 
fluid dynamical baseline for w 2 as a function of collision energy, centrality and system 
size, measurements of this suppression could thus be used to constrain the QGP 
shear viscosity from experimental data. Other observables (e.g. slopes of 7>r-spectra) 
also depend on 77/s, but i>2 appears to show the strongest sensitivity, so it has been 
studied most extensively, and we will focus on it. We will concentrate on viscous 
hydrodynamics, leaving a discussion of more microscopic approaches (e.g. [11, 12]) for 
another time. 

The first results for viscous v-i suppression published by the different groups 
seemed to show large discrepancies, ranging from 20% to 70% for 'minimal viscosity' 
77/s = 1/47T. A detailed code verification process [13], carried out within the TECHQM 
collaboration [15], eliminated the possibility that these differences were caused by 
numerical error. Systematic studies during the last few months revealed their physical 
origins and are summarized in this talk. Specifically, we will discuss the effects of 
system size, equation of state (EOS), and different versions of the Israel-Stewart 
equations used to evolve the viscous terms in the energy-momentum tensor, on 
viscous V2 suppression. A combination of these effects has been shown to resolve the 
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apparent discrepancies reported by different groups [6]. Further, we will discuss effects 
from different freeze-out procedures [16], different initializations (Glauber model vs. 
Color Glass Condensate approach) [7], bulk viscosity [17], and numerical viscosity [6]. 
Quantifying these effects is important for understanding the uncertainties in extracting 
the shear viscosity to entropy ratio from experimental data. 

2. Viscous hydrodynamics in 2+1 dimensions 

In this section, we briefly review the causal viscous equations solved by the VISH2+1 
code (Viscous Isreal-Stewart Hydrodynamics in 2+1 space-time dimensions) developed 
at OSU, used to simulate QGP and hadronic matter expansion with exact longitudinal 
boost invariance but arbitrary dynamics in the transverse directions [18]. For 
simplicity, we assume zero net baryon number and heat conductivity. VISH2+1 solves 
the equations for energy momentum conservation d m T mn = 0, with 

T mn = eu m u n _ ^ + n ) A mn + ^ran ^ ^ran = ^mn^r^n (!) 

together with kinetic equations for the viscous shear pressure tensor ir mn and the bulk 
pressure IT: 

A mr A ns Dir rs = - — (7T mn -2 W mn ) - l -Ti mn ^-d k ( —u k ) , (2) 
t v 2 tv \r)T ) 

DU = --(H + (9) \u C —d k (^u k ) . (3) 
T n 2 7- n \Q1 ) 

Here D = u m d m is the time derivative in the local comoving frame, V m = A ml di is the 
spatial gradient in that frame, and a mn = V<" 1 w™> = \{V m u n +V n u rn ) - \A mn 6 (with 
the scalar expansion rate 9 = d k u k = V k u k ) is the symmetric and traceless velocity 
shear tensor. d m denotes the components of the covariant derivative in the curvilinear 
coordinates (T 7 x,y,ri s ) [18]. 

For systems with conformal symmetry the last terms in Eqs. (2,3) can be written 
in different forms as discussed in [6]. While these forms are equivalent for conformal 
systems, they differ in principle for systems with an EOS that breaks conformal 
invariance, e.g. through a phase transition. These differences turn out to be negligible 
in practice [6]. Without the last terms, Eqs. (2,3) are known as simplified IS equations 
whereas their complete versions will be called full IS equations. In spite of this 
name, these equations still do not include all possible second-order terms in a gradient 
expansion around the ideal fluid (locally thermalized) limit. For a conformal theory 
with vanishing chemical potentials 4 other terms can be added to the right hand side of 
Eq. (2), with additional coefficients Ai, A 2 , A 3 , and n (k = in Minkowski space) [19]. 
These terms include couplings to the vorticity tensor [19] which turns out to be small 
in heavy-ion collisions if the initial longitudinal velocity profile is boost invariant [1]. 
Even more terms arise in a kinetic theory derivation that does not assume conformal 
symmetry [20], including terms that couple Eqs. (2) and (3). Their coefficients can 
be obtained from kinetic theory [20, 21] at weak coupling or from the AdS/CFT 
correspondence at strong coupling [19]. Still another approach was developed by 
Ottinger and Grmela (see references in [3]); when translated into Israel-Stewart form, 
it also makes specific predictions for the coefficients of these second-order terms [16]. 
The so far accumulated numerical evidence [6, 7, 13] suggests that, except for the last 
terms in Eqs. (2) and (3), all other second order terms are unimportant in practice, 
but a systematic study that confirms this beyond doubt remains outstanding. 
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The shear viscosity 77, bulk viscosity £ and the corresponding relaxation times 
TV and Tn are free parameters in VISH2+1. For most of the published numerical 
simulations of (2+l)-d viscous hydrodynamics, the shear viscosity has been set to 
the KSS minimal value r//s = 1/Air ~ 0.08, and the bulk viscosity was set to zero, 
C/s = 0. Bulk viscous effects on the evolution of elliptic flow in non-central collisions 
were recently studied by us [17] and are here reported for the first time (see Sec. 5). 
For the relaxation times tv and rn most authors used a constant multiple of r\j sT . 
Physical observables turn out to be largely insensitive to the value of tv if (and only 
if!) the last terms are included in Eqs. (2) and (3), i.e. the full (and not the simplified!) 
Israel-Stewart equations are used [6] (see Sec. 3 below). 

The equation of state (EOS), initial and final conditions are additional inputs for 
both ideal and viscous hydrodynamics. Details can be found in the original literature 
[3, 4, 6, 7] and will not be explained here. Two different EOS will be used in these 
proceedings: SM-EOS Q is a slightly smoothed version [4] of EOS Q [22] which describes 
a non-interacting massless QGP phase matched to a chemically equilibrated massive 
hadron resonance gas (HRG) through a Maxwell construction. EOS L matches the 
HRG EOS below T c smoothly with the lattice QCD EOS above T c [4]. Our EOS L is 
close to but not identical with the quasiparticle EOS used by Romatschke [1, 7]. 

3. Viscous V2 suppression: effects from system size, EOS and different 
versions of Israel-Stewart equations 

In this section we will briefly discuss the different manifestations of shear viscosity 
when one varies system size and EOS and uses different versions of the I-S equations 
[6]. As mentioned in the Introduction, this analysis resolves the initially puzzling 
differences between the results published by different groups. 
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Figure 1. (Color online) Differential elliptic flow V2{pt) f° r directly emitted 
pions (i.e. without resonance decay contributions), comparing results for different 
collisions systems and EOSs, for ideal and viscous fluid dynamics, with parameters 
as indicated. 

Figure 1 shows the differential elliptic flow V2(pt) for directly emitted pions from 
ideal and viscous hydrodynamics. Panels (a) and (b) compare two systems of different 
size (Cu+Cu and Au+Au, both at b = 7fm), using the same equation of state (SM- 
EOS Q), I-S equations (simplified) and other free inputs. Although both systems have 
similar initial eccentricities, the smaller Cu+Cu system shows a much larger viscous 
V2 suppression (by almost 70% below the ideal fluid value at pr = 2GeV/c [2, 4]) 
than observed in the larger Au+Au system where the suppression is almost a factor 
two smaller. Panels (b) and (c) compare the same Au+Au system at b = 7fm for 
two different EOS and different I-S equations. Changing the EOS from SM-EOS Q to 
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EOS L reduces the viscous suppression of elliptic flow by another quarter (from ~ 40% 
to ~ 30% below the ideal fluid value at px = 2GeV/c). Replacing the simplified I-S 
equations used in [2, 4] by the full I-S equations used in [1] further reduces the v 2 
suppression from 30% to 25% below the ideal hydrodynamics value at pt = 2GeV/c. 
This final result is consistent with [1]. Although for EOS L the additional term in the 
full I-S equations only results in a 5-10% difference in v 2 suppression, its effect is much 
larger for more rapidly expanding systems, such as Cu+Cu collisions driven by a stiff 
conformal EOS e = 3p [6]. In such systems the simplified I-S equations lead to a strong 
sensitivity of physical observables on the value of the microscopic relaxation time T n , 
whereas inclusion of the last term in the full I-S equations (2) appears to minimize this 
sensitivity in all situations [6]. (It is also needed to preserve the conformal symmetry 
in conformal systems [19].) 

4. Dynamical freeze-out and effects from late hadronic viscosity 

In ideal hydrodynamics, one usually imposes "sudden freeze-out", i.e. a sudden 
transition from thermalized fluid to free-streaming particles on a hypersurface S(x) 
of constant temperature or energy density [22]. The same algorithm has also been 
used in most of the existing viscous hydrodynamic calculations [1, 2, 4, 5, 6, 7]. 
Since viscous hydrodynamics is based on an expansion in small deviations from local 
equilibrium, its validity requires the microscopic relaxation time to be much smaller 
than the macroscopic inverse expansion rate, r ve id-u <C 1. This condition (whose 
long history is discussed in Ref. [23] where it is also applied to ideal hydrodynamics) 
provides a natural criterium for a dynamical freeze-out algorithm. Dusling and Teaney 
[3] implemented it into their viscous hydrodynamics. They find that in this case the 
viscous v 2 suppression is dominated by non-equilibrium corrections to the local thermal 
distribution function along the freeze-out surface [3]. This is not a collective effect 
arising from anisotropics of the flow velocity profile, but a reflection of non-equilibrium 
momentum anisotropics in the local fluid rest frame. In contrast, for isothermal freeze- 
out we find [4] that the viscous v 2 suppression is dominated by the viscous reduction 
of the collective flow anisotropy, while local rest frame momentum anisotropics play 
a much smaller role. This comparison shows that a careful treatment of the hadronic 
decoupling process will be required for the quantitative extraction oir)/s from elliptic 
flow data. Dusling also found that dynamical freeze-out can increase the slope of the 
multiplicity dependence of the eccentricity-scaled elliptic flow v 2 /e [16]. This is an 
improvement over the scaling behaviour found in [6] for viscous hydrodynamics with 
constant rj/s and isothermal freeze-out which features a slope that is too small. 

There are other reasons why a proper kinetic treatment of the late hadronic 
phase is important. By matching a realistic hadron rescattering cascade to an ideal 
fluid description of the QGP and hadronization stages, Hirano et al. [24] showed that 
the HRG phase is highly viscous and strongly suppresses any buildup of elliptic flow 
during the hadronic stage. This is consistent with a recent analysis by Demir and 
Bass [25] who found large shear viscosities for their hadronic UrQMD cascade even 
close to T c (between 5 and 10 times above the KSS bound). For a correct description 
of the beam energy and centrality dependence of v 2 , which are crucially affected by 
the changing relative weight of QGP and HRG dynamics in building elliptic flow 
(in central collisions or at higher energies the system spends more time in the QGP 
phase than in peripheral collisions or at low energies), a realistic kinetic simulation 
of the hadronic phase and its freeze-out thus appears to be indispensable. It will 



Extracting QGP viscosity from RHIC data 



5 



also solve the problem of using an incorrect chemical composition in the hadronic 
phase that plagues present implementations of viscous hydrodynamics: the empirical 
finding at RHIC that chemical equilibrium is broken close to T c and does not (as so far 
assumed in all viscous codes) persist to the point of kinetic freeze-out has important 
consequences for elliptic flow. As shown in Refs. [26, 27] within the framework of ideal 
fluid dynamics, the distribution of the total momentum anisotropy among the various 
hadronic species depends strongly on the chemical composition of the hadronic fireball 
at freeze-out, with almost 25% larger pion elliptic flow when chemical equilibrium 
is broken at T c ss 165 MeV than for the case where pions are allowed to remain 
in chemical equilibrium down to Td oc ~ 100 MeV. It is obvious that such a large 
effect must be correctly implemented in viscous hydrodynamics before a quantitative 
extraction of rj/s for the QGP can be attempted. 

5. Bulk viscosity 
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Figure 2. (Color online) Left: The bulk viscosity to entropy ratio as a function 
of temperature, as used in these proceedings (see text for discussion). Right: 
Differential elliptic flow V2(pt) for directly emitted pions calculated with ideal 
hydrodynamics and with viscous hydrodynamics. The viscous calculations assume 
minimal shear viscosity and three different normalizations of the bulk viscosity 
shown in the left panel, as noted in the legend. 

Early viscous hydrodynamic simulations ignored bulk viscosity for simplicity. 
Although it vanishes for a conformal fluid or massless QGP on the classical level, 
quantum effects break the conformal symmetry of QCD and generate a nonzero 
bulk viscosity even in the massless QGP phase, as recently measured using lattice 
QCD [28]. General arguments predict that the specific bulk viscosity (/s exhibits 
a peak near T c [29, 30] and then decreases again on the hadronic side. However, 
different theoretical approaches show dramatically different peak values: the one from 
AdS/CFT predictions [31] is more than a order of magnitude smaller than the values 
extracted from lattice QCD data [28] (see, however, Ref. [32] for a critical discussion 
of the lattice QCD approach). Considering this theoretical uncertainty, we treat the 
bulk viscosity as an essentially free input, implementing only the feature that it peaks 
around T c . To obtain an explicit expression for (£/s)(T) we connect the minimum 
AdS/CFT result from Ref. [33], evaluated with lattice QCD data for the speed of sound 
c^(T), through a Gaussian function peaked at T c with a zero value in the hadronic 
phase (lower brown line (C = 1) in the left panel of Figure 2). To simulate effects 
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from larger bulk viscosity, we multiply the entire function by a factor C > 1 (C = 2 
for the upper line in the left panel of Figure 2). 

Since the fireball expands, the Navier-Stokes limit of the bulk viscous pressure 
II = — (8 is negative. This reduces the thermal pressure (effectively softening the EOS 
near T c ), decelerates longitudinal expansion and suppresses the buildup of radial flow. 
As a result, the hadron p-r-spectra become steeper. By construction, the bulk viscous 
effects arc strongest around T c , but since the transverse edge of the fireball already 
hadronizes early when the longitudinal expansion rate is high, this does not imply 
that bulk viscous effects are necessarily small at early times. Indeed, we see in the 
right panel of Fig. 2 that elliptic flow is significantly reduced by bulk viscosity even 
though a large fraction of the momentum anisotropy is generated well before most of 
the matter reaches T c . For minimal shear viscosity r\js = 1/47T, we find that bulk 
viscosity can increase the viscous suppression of pion elliptic flow at pr = 0.5 GeV by 
an additional 25% (for C = 1) to 50% (for C = 2) relative to the case C, = 0. 

Clearly, bulk viscous effects must be accounted for when trying to extract n/s 
from elliptic flow data. At this moment it is not clear which observables can be used 
for a clean separation of bulk and shear viscous effects. One should explore whether 
the centrality dependence of the local hydrodynamic expansion rate impacts radial 
and elliptic flow in distinct ways that allow for such a separation. 

6. Glauber model vs. CGC initialization 

We now come to what may turn out to be a serious road block for precision 
measurements of the QGP shear viscosity: our insufficient knowledge of the initial 
source eccentricity e. It has now been known for while that the Color Glass Condensate 
(CGC) model, implemented in the initial entropy or energy density profile via the 
KLN parametrization (see [24, 34] for references), leads to ~ 30% larger initial source 
eccentricities than the popular Glauber model. Ideal fluid dynamics transforms this 
larger source eccentricity into ~ 30% larger elliptic flow. Since the extraction of n/s is 
based on the viscous suppression of v 2 , obtained by comparing the measured elliptic 
flow with an ideal (inviscid) fluid dynamical benchmark calculation, a 30% uncertainty 
in this benchmark can translate into a 100% uncertainty in the extracted value for 
n/s. This was recently shown by Luzum and Romatschke (see Fig. 8 in [7]). This 
uncertainty trumps most of the other uncertainties discussed above. Worse, since the 
initial source eccentricity depends on details of the shape of the KLN profile near 
the edge of the distribution where the gluon saturation momentum scale Q s becomes 
small and the CGC model reaches its limit of applicability, there is little hope that we 
can eliminate it theoretically from first principles. 

Based on their analysis of charged hadron elliptic flow data from the STAR 
experiment, allowing for a 20% systematic uncertainty of these data, the authors of 
[7] found an allowed range < n/s < 0.1 for Glauber and 0.08 < n/s < 0.2 for CGC 
initial conditions. Since their analysis did not include a comprehensive investigation 
of effects caused by permissible variations of the EOS near T c , by bulk viscosity, or by 
late hadronic viscosity and non-equilibrium chemical composition at freeze-out (see 
Sees. 3-5 above), one should add a significant additional uncertainty band to these 
ranges. Furthermore, correcting the experimental data for event-by-cvent fluctuations 
in the initial source eccentricity [35] may bring down the measured V2 values even 
below the range considered in [7]. Still, we agree with Luzum and Romatschke that, 
even when adding all the above effects in magnitude (ignoring the fact that several 
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of them clearly have opposite signs), viscous hydrodynamics with rj/s > 5 x (l/47r) 
would suppress the elliptic flow too much to be incompatible with experiment. 

7. Numerical viscosity 

To study the effects from shear and bulk viscosity one must ensure that numerical 
viscosity is under control and sufficiently small. Simply speaking, numerical viscosity 
comes from the discretization of the hydrodynamic equations for numerical calculation. 
It causes entropy production even in ideal hydrodynamics without shocks and can 
never be fully avoided. To minimize numerical viscosity, the flux-corrected transport 
algorithm SHASTA [36] employed by VISH2+1 (and by its ideal fluid ancestor AZHYDRD 
[14]) implements an "antidiffusion step" involving a parameter A called "antidiffusion 
constant" [36] . For a given grid spacing, numerical viscosity is maximized by setting 
A = 0. In standard situations, the default value A = g minimizes numerical viscosity 
effects [36]. With A = | and typical grid spacing Ax = Ay = 0.1 fm, At = 0.04 fm/c, 
AZHYDRO generates only 0.3% additional entropy in central Au+Au collisions. This is 
negligible when compared with the O(10%) entropy production by VISH2+1 for a fluid 
with real shear viscosity n/s — 1/Aw. 

By increasing the grid spacing in AZHYDRO and/or changing A, we can explore 
the effects of numerical viscosity on radial and elliptic flow. We find that numerical 
viscosity has little effect on the development of radial flow but reduces V2 in very much 
the same way as does real shear viscosity. Since we gauge the effects of rj/s on vi by 
comparing results from VISH2+1 for rj/s ^ to those for r\j s = 0, we should explore 
how much in the latter case is already suppressed by numerical viscosity. We can 
do this by setting r]/s = and reducing the grid spacing until v 2 stops changing (i.e. 
until we have completely removed all numerical viscosity effects on v%). In this way 
we have ascertained that for our standard grid spacing numerical viscosity suppresses 
the differential elliptic flow v 2 (pt) by less than 2%. 

8. Summary and outlook 

While the elliptic flow u 2 generated in non-central heavy-ion collisions is very sensitive 
to the shear viscosity to entropy ratio n/s of the QGP, it is also significantly affected 
by (i) details of the initialization of the hydrodynamic evolution, (ii) bulk viscosity 
and sound speed near the quark-hadron phase transition, and (iii) the chemical 
composition and non-equilibrium kinetics during the late hadronic stage. Not all of 
these effects are presently fully under control. Recent attempts to extract the specific 
shear viscosity rj/s phenomenologically, by comparing experimental elliptic flow data 
with viscous hydrodynamics, have established a robust upper limit 

- <5xp (4) 

S QGP 47T 

close to the conjectured KSS bound [10], but further progress requires elimination 
of the above systematic uncertainties. Since some of these influence the build-up 
of elliptic flow in opposite directions, it is quite conceivable that the QGP specific 
viscosity is in fact much closer to the KSS bound (77/ s)| KSS = 1/4-7T than suggested by 
the upper limit (4). Ongoing improvements on the theory side should help to reduce 
or eliminate most of the above uncertainties, bringing us closer to a quantitative 
extraction of rj/s for the quark-gluon plasma. The single largest uncertainty, however, 
is caused by our poor knowledge of the initial source eccentricity which varies by about 
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30% between models. As shown in [7], this translates into an 0(100%) uncertainty 
for rj/s. It seems unlikely that theory can help to eliminate this uncertainty from first 
principles. It thus appears crucial to develop experimental techniques that may help 
us to pin down the initial source eccentricity phenomenologically, with quantitative 
precision at the percent level. 
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